Fusariose: QTLs Description - Part2

7/ IC on the map (Fig)
8/ Evolution of LOD scores with time?
- Functions
- QTL 1B
- QTL 5A
9/ Link mith physical reference
- 1B
- 5A
- 4B (height)
- 7B (precocity)
- summary
10/ Comparaison avec les QTLs de la biblio
- Using SSR
- Position of Snn1 on map
11/ Relationship QTL presence / field characteristic
12/ Epistatic effect?
- Example
- All trait

This file follows the first part of the QTL analysis (Description_of_QTLs-Part1). Load the results + functions of the first part:

load("/Users/yan/Dropbox/Publi_Fusariose/ANALYSIS_REPRO/DATA/QTL_R_environment.R")

Charge some libraries that will be useful

library(RColorBrewer)
library(xtable)
library(plotly)
library(FactoMineR)

7/ IC on the map (Fig)

We can also visualize the IC on the genetic map. Let’s visualize important variable ONLY = AUDPC and last note for BLUPs only.

par(mfrow=c(1,1))
my_var_to_show=PUBLI
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)

I save it as a png image for publication purpose.

png("~/Dropbox/Publi_Fusariose/FIGURE/Fig_IC_QTL_on_map.png")
show_IC_QTL_on_map(data, my_var_to_show , c("1B", "5A"), my_QTL_thres , 1.5)
dev.off()

## quartz_off_screen 
##                 2

8/ Evolution of LOD scores with time?

We are going to study qtl 1B and 5A only. Let’s have a look to the maximum LOD score of these QTL at each date of observation. For each experiment, for each variable?

Functions

Fonction qui regroupe chaque LOD max dans un tableau. On ne garde que les BLUPs de toutes les variables à toutes les dates.

find_LOD_max=function( chromo, my_min, my_max   ){
  bilan=data.frame(matrix(0, 200  ,8))
  num=0
  vec_names=c("NOT", "PEPI", "NEPIL", "PEPIL")
  for(expe in list(CAP13, GRI11, GRI13, GRI15, LEC14, BEQ11, BAR14)){
    
    num_carac=0
    for(carac in list(NOT, PEPI, NEPIL, PEPIL)){
      num_carac=num_carac+1
      my_carac=intersect(BLUP, intersect(expe, carac ))
      if( expe==BEQ11){ my_carac=intersect(expe, carac ) }
      if( length(my_carac)==0) { next }
      num=num+1
      my_vec=c()
      for(i in c(1:length(my_carac))){
        b=(data$LOD[ which(data$LG==chromo & data$variable==my_carac[i] & data$Distance<my_max & data$Distance>my_min ) ])
        a=max( b , na.rm=T)
        my_vec=c(my_vec, a)
      }
      
      my_replace=c( unique(gsub("_.*","",my_carac)) , vec_names[num_carac], my_vec )
      bilan[num, c(1:length(my_replace)) ] = my_replace 
    }
  }
  bilan=bilan[bilan[,1]!=0 , ]
  bilan[bilan==0]=NA
  return(bilan)
}

Fonction pour visualiser ce tableau.

show_max_LOD=function(bilan){
  
  # Préparation des axes des X.
  x1=c("300", "400", "550", "AUDPC")
  x2=c("300", "350", "450", "500", "550", "AUDPC")
  x3=c("350", "450", "550", "AUDPC")
  x4=c("350", "550", "AUDPC")
  x5=c("300", "400", "550", "AUDPC")
  x6=c("300", "350", "450", "500", "AUDPC")
  x7=c("300", "400", "550", "600", "AUDPC")
  my_x_axis=list(x1,x2,x3,x4,x5,x6,x7)
  
  par(mfrow=c(2,4) , mar=c(3,3,1,1))
  liste_var=c("NOT", "PEPI","NEPIL", "PEPIL")
  num=0
  for(expe in list("CAP13", "GRI11", "GRI13", "GRI15", "LEC14", "BEQ11", "BAR14")){
    num=num+1
    AA=bilan[ which(bilan[,1]==expe),]
    plot(seq(1,6) , AA[ 1, c(3:ncol(AA))] , type="l", ylim=c(0,10) , col=my_colors[match( AA[ 1, 2] , liste_var )], lwd=2, xlab=expe, ylab="", xaxt="n")
    axis(1, at=seq( 1 , length(my_x_axis[[num]])), my_x_axis[[num]], cex=0.7)
    abline(h=my_QTL_thres, col="grey")
    abline(v=c(1:6), col="grey", lwd=0.3)
    for(i in c(2:nrow(AA))){
     points(seq(1,6) , AA[ i, c(3:ncol(AA))] , type="l" , col=my_colors[match( AA[ i, 2] , liste_var )], lwd=2)  
    }
   text(x=1.8, y=9, expe, col="orange" )
  }
  plot(1,1,col="transparent", xaxt="n", yaxt="n", bty="n", xlab="", ylab="")
  legend( "center", legend = c("NOT", "PEPI","NEPIL", "PEPIL") , bty="n", pch=20 , col=my_colors[1:4], pt.cex=4, cex=1.6 )
}

QTL 1B

Calcul des LOD max à chaque date:

max_lod_1B=find_LOD_max( "1B", 0, 9 )

Show it:

show_max_LOD(max_lod_1B)

Donc ces graphiques représentent les valeurs des LODs max de chaque QTL, i.e. pour chaque variable de chaque expe pour le chromosome 1B.

QTL 5A

Calcul des LOD max à chaque date:

max_lod_5A=find_LOD_max( "5A", 100, 300   )

Show it:

show_max_LOD(max_lod_5A)

9/ Link mith physical reference

OK so we have a QTL with a confidence interval. It could be great to see how many genes are present in this IC, and the physical size it represents.

1B

Summary of the 1B IC of the QTL:

# initialise the table
bil_phy=data.frame(matrix(0,0,9))
colnames(bil_phy)=c("chromo","start IC (cM)", "end IC (cM)", "len IC (cM)", "# SNP", "start IC (Mb)", "end IC (Mb)", "size IC (Mb)", "# genes")

# 1B
bil_phy[1,]=find_IC_infos(1.2,14.1,"1B",c())

chromo	start IC (cM)	end IC (cM)	len IC (cM)	# SNP	start IC (Mb)	end IC (Mb)	size IC (Mb)	# genes
1B	1.2	14.1	12.5	58	1527986	6092666	4564680	0

This IC has 103 genes. Let’s check the colinearity between genetic and physical positions in this area.

#1B
par(mfrow=c(1,2), mar=c(5,5,3,3))
show_genes_of_IC(0, 300, "1B", 0.005, 1, c())
show_genes_of_IC(0, 14, "1B", 0.5, 4, c())

5A

Summary of the 5A IC of the QTL:

#5A
bil_phy[2,]=find_IC_infos(235.6 ,305,"5A","Traes_5AL_99116D900@696")

chromo	start IC (cM)	end IC (cM)	len IC (cM)	# SNP	start IC (Mb)	end IC (Mb)	size IC (Mb)	# genes
5A	235.6	305	68.7	86	118853998	143979436	25125438	0

#5A
par(mfrow=c(1,2), mar=c(5,5,3,3))
show_genes_of_IC(0, 400, "5A", 0.005, 1, c())
show_genes_of_IC(235.6 ,305, "5A", 0.1, 4,"Traes_5AL_99116D900@696")

4B (height)

Summary of the 4B IC of the QTL:

# 4B
bil_phy[3,]=find_IC_infos(47,53,"4B","")

## Warning in min(IC_map$Distance): aucun argument trouvé pour min ; Inf est
## renvoyé

## Warning in max(IC_map$Distance): aucun argument pour max ; -Inf est renvoyé

## Warning in min(IC_map$Posi_physique): aucun argument trouvé pour min ; Inf
## est renvoyé

## Warning in max(IC_map$Posi_physique): aucun argument pour max ; -Inf est
## renvoyé

chromo	start IC (cM)	end IC (cM)	len IC (cM)	# SNP	start IC (Mb)	end IC (Mb)	size IC (Mb)	# genes
4B	47	53	-Inf	0	Inf	-Inf	-Inf	0

#4B
par(mfrow=c(1,2), mar=c(5,5,3,3))
#show_genes_of_IC(0, 200, "4B", 0.005, 1)
#show_genes_of_IC(47, 53, "4B", 0.1, 4)

7B (precocity)

Summary of the 7B IC of the QTL:

# 7B
bil_phy[4,]=find_IC_infos(5,35,"7B","")

## Warning in min(IC_map$Distance): aucun argument trouvé pour min ; Inf est
## renvoyé

## Warning in max(IC_map$Distance): aucun argument pour max ; -Inf est renvoyé

## Warning in min(IC_map$Posi_physique): aucun argument trouvé pour min ; Inf
## est renvoyé

## Warning in max(IC_map$Posi_physique): aucun argument pour max ; -Inf est
## renvoyé

chromo	start IC (cM)	end IC (cM)	len IC (cM)	# SNP	start IC (Mb)	end IC (Mb)	size IC (Mb)	# genes
7B	5	35	-Inf	0	Inf	-Inf	-Inf	0

#7B
#par(mfrow=c(1,2), mar=c(5,5,3,3))
#show_genes_of_IC(0, 200, "7B", 0.005, 1)
#show_genes_of_IC(5, 35, "7B", 0.1, 4)

summary

chromo	start IC (cM)	end IC (cM)	len IC (cM)	# SNP	start IC (Mb)	end IC (Mb)	size IC (Mb)
1B	1.2	14.1	12.5	58	1527986	6092666	4564680
5A	235.6	305	68.7	86	118853998	143979436	25125438
4B	47	53	-Inf	0	Inf	-Inf	-Inf
7B	5	35	-Inf	0	Inf	-Inf	-Inf

10/ Comparaison avec les QTLs de la biblio

Using SSR

Let’s have a look to all the 9 SSR markers I have in my map:
Is it normal I have only 9?
This is not gonnna help a lot…
- 1B: pas de SSR
- 5A: 2 SSR, mais très très loin de l’IC.

–> Il faut aller voir si je peux trouver des arqueurs de la littérature qui ont été passé sur Dic2 x Silur, et les replacer sur la carte comme pour la publi mosaiquE…

SSR=unique(data[ grep("gwm" , data$marqueur) , c(1:3) ])
nrow(SSR)

## [1] 9

LG	marqueur	Distance
1A	Xgwm164	68.40
3B	Xgwm376	92.30
3B	Xgwm77	94.20
4B	Xgwm495	94.00
5A	Xgwm205	41.00
5A	Xgwm415	72.50
5B	Xgwm544	55.30
7A	Xgwm276	181.50
7B	Xgwm400	38.50

Position of Snn1 on map

Le gène Snn1 est positionné sur le 1B et apporte une résistance forte à la septoriose, un autre champignon. On peut trouver la séquence de ce gène ici:
https://www.ncbi.nlm.nih.gov/nuccore/858936371/ Publié par: http://advances.sciencemag.org/content/2/10/e1600822.full

Blast sur la référence blé tendre release 28

makeblastdb -in sequence_snn1.fasta -dbtype nucl -out database

Blast à 100% sur ce contig: Traes_1BS_C4B727779

Ce contig n’est pas dans notre carte génétique, donc on a pas de SNP dessus. Par contre, en terme de position physique, il se situe à 880274 Mb sur le 1BS. Ca correspondrait donc à être entre 1.19 et 1.59 cM sur notre carte.

1B Traes_1BS_7BC6EFB55@1326 1.19999999999999 1B 502322 1B Traes_1BS_81F7DAE3E@432 1.59999999999999 1B 2538668

Nous notre QTL a tendance a etre un poil plus loin, mais vraiment vraiment proche…

11/ Relationship QTL presence / field characteristic

The intensity of the Fusarium infection was not the same from one experiment to another. Let’s check the relationship between disease intensity and strength of the QTL detected?

12/ Epistatic effect?

OK so we have 2 main QTLs. Do we have a epistatic interaction betwwen them??

Example

Let’s start on an example: BAR14_PEPIL300_blup. 2 QTLs: markers Cluster_16778|Contig1|original@1840 et Cluster_9940|Contig1|complementarySeq@104

# calculate inter
bil_inter=data.frame(matrix(0,0,9))
colnames(bil_inter)=c("trait","marker-1B", "marker 5A", "mean A-A", "mean A-B", "mean B-A", "mean B-B", "R2 tot", "pval inter")
bil_inter[1,]=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")

trait	marker-1B	marker 5A	mean A-A	mean A-B	mean B-A	mean B-B	R2 tot	pval inter
BAR14_PEPIL300_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	-0.2	3.32	-0.91	-0.65	0.281	0.0001369779

And show this interaction!

boxplot_two_QTL_interactive("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup")

#boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , "BAR14_PEPIL300_blup", ylim=c(-4,7))

All trait

Il faut commencer par trouver tous les traits qui ont un QTL signif sur le 1B et un signif sur le 5A.

var_with_2QTLs=c()
for(i in list(summary_QTL_NOT, summary_QTL_PEPI, summary_QTL_NEPIL, summary_QTL_PEPIL)){
  a=i[which(i$chromo%in%c("1B","5A")) , ]
  b=table(a$carac)
  var=names(b[b==2])
  var_with_2QTLs=c(var_with_2QTLs, var)
}
var_with_2QTLs=var_with_2QTLs[grep("blup",var_with_2QTLs)]

On a 8 variables qui ont les 2 QTLs significatifs. Voila de quelles variables il s’agit:

var_with_2QTLs
BAR14_PEPI350_blup
BAR14_PEPI400_blup
BAR14_PEPI600_blup
BAR14_PEPIAUDPC_blup
GRI13_PEPI550_blup
BAR14_NEPIL300_blup
BAR14_PEPIL300_blup

Calculons les p-values de l’effet d’interaction pour tous ces cas:

for(i in var_with_2QTLs[-which(var_with_2QTLs=="BAR14_PEPIL300_blup")]){
  vv=analyse_inter("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , i)
  bil_inter=rbind(bil_inter , vv)
}

trait	marker-1B	marker 5A	mean A-A	mean A-B	mean B-A	mean B-B	R2 tot	pval inter
BAR14_PEPIL300_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	-0.2	3.32	-0.91	-0.65	0.281	0.0001369779
BAR14_PEPI350_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	0.38	21.58	-7.91	-3.92	0.355	0.0004079017
BAR14_PEPI400_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	-0.08	20.57	-7.65	-4.53	0.244	0.0039045779
BAR14_PEPI600_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	2.17	3.69	-4.75	1.02	0.318	0.0163055323
BAR14_PEPIAUDPC_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	896.84	9825.26	-4380.44	-2080.69	0.394	0.002065465
GRI13_PEPI550_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	-1.25	0.94	0.54	1.12	0.023	0.4583633636
BAR14_NEPIL300_blup	Cluster_16778\|Contig1\|original@1840	Cluster_9940\|Contig1\|complementarySeq@104	-1.33	18.18	-4.85	-3.36	0.259	0.0002526347

On a donc un effet d’épistasie qui est significatif à tous les coups (quasiment).

Représentons ces épistasies avec un boxplot pour chacune des variables:

# widget to choose data
#selectInput("select", label = h3("Select a variable"), choices = var_with_2QTLs,  selected = 1)

# interactibe boxplot
#renderPlotly({
#  boxplot_two_QTL("Cluster_16778|Contig1|original@1840", "Cluster_9940|Contig1|complementarySeq@104" , select)
#})

Dans la majorité des cas, on observe que les individus avec 1 seul QTL ou avec 2 QTLs sont similaires. Par contre, les individus avec aucun des 2 sont très différents et très sensibles. Donc en gros, si j’ai au moins un des 2 gènes tout va bien, sinon je suis sensible. Donc on est pas dans un modèle additif à priori, les 2 gènes sont liés.